/*
 * Moseycode
 * Copyright (C) 2008  Tom Gibara
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#include <limits.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <jni.h>
#include <android/log.h>



/* Some quick and dirty macros to make logging convenient */
#define LOG_TAG "mosey.cpp"
#define LOGE(...) ((void)__android_log_print(ANDROID_LOG_ERROR, LOG_TAG, __VA_ARGS__))
#define LOGI(...) ((void)__android_log_print(ANDROID_LOG_INFO, LOG_TAG, __VA_ARGS__))

/* Control debug */

// defining this causes the PLOT array to accumulate an informative greyscale image
#define MOSEY_DEBUG

/*
 * SHARED CONSTANTS
 *
 * These parameters constrain the processing and are fixed at compile time.
 */

// the maximum frame sizes that can be processed
// nothing prevents these maximums from being increased
// except that memory allocations are statically based on MAXSIZE
static const int MAXWIDTH = 256;
static const int MAXHEIGHT = 256;
//TODO decouple definition of MAXSIZE from MAXWIDTH and MAXHEIGHT
static const int MAXSIZE = MAXWIDTH * MAXHEIGHT;

/*
 * SHARED SETTINGS
 *
 * These parameters are set during initialization and control the characteristics
 * of the frame processing. The values can be reinitialized to modify the parameters.
 */

//determines extent of averaging top adapt to variations in brightness
//the intensity average is taken over 2^ORDER pixels horizontally and vertically
static int ORDER;
// the width in pixels of the frames to be processed
static int WIDTH;
// the height in pixels of the frames to be processed
static int HEIGHT;

/*
 * All values under a DERIVED haeding are derived from the shared settings above.
 * They are recomputed on each initialization. All values under a CONSTANT heading
 * are independent of the shared settings.
 */

/*
 * SHARED DERIVED
 *
 * Computed values relevant to all functions.
 */

//the total number of pixels: WIDTH * HEIGHT
static int SIZE;

/*
 * AVERAGING CONSTANTS
 *
 * Fixed parameters for the average function.
 */

//TODO make these parameterizable
static const int BIAS = 7;
static const int BIAS_ORDER = 3;

//TODO initialization not necessary?
static int AVERAGES[MAXSIZE] = { 0 };

/*
 * AVERAGING DERIVED
 *
 * Precomputed values for the average function.
 */

//dimensions of filter
static int HORDER;
static int VORDER;

//horizontal and vertical sizes
static int HSIZE;
static int VSIZE;
static int SWIDTH;
static int SHEIGHT;

//offsets for locating oldest summand
static int HOFFSET;
static int VOFFSET;

//half horizontal and vertical sizes
static int HHSIZE;
static int HVSIZE;
static int HSWIDTH;
static int HSHEIGHT;

//the amount we ultimately need to shift by in order to recover the avg.
static int SHIFT;

/*
 * DETECTOR CONSTANTS
 *
 * Fixed parameters for the detect function.
 */

// Whether the neutral 'color' is true (light) or false (dark).
static const bool NEUTRAL = true;

//enumeration of directions
static const char E = 0;
static const char S = 1;
static const char W = 2;
static const char N = 3;

//value assigned when not an edge pixel
//TODO bad assumes unsigned chars
static const char BAD_PLACE = 255;

//the maximum length of a path before the algorithm gives up
//it's computed as the length of the frame perimeter
//(this is the largest convex perimeter possible)
static const int MAXLEN = 4 * (MAXWIDTH + MAXHEIGHT);
//the minimum length of a path that will be considered for shape analysis
static const int MINLEN = 6;

//an array which records which 'north' path movements have been made
static bool USED[MAXSIZE];

//the number of path steps desirable between points on a path
//this is specified as a power of two for fast division
//and the result of the division is bounded by the specified the max/min
static const int PTSPACE_ORDER = 3;
static const int MAXPOINTS = 12;
static const int MINPOINTS = 4;
//TODO should this be based on width&height?
static const int NONSIG = 20;

//accumulates directions for the current path
static char PATH[MAXLEN];
//records list of integer coordinates along the current path
static int POINTS[2*MAXPOINTS];

/*
 * DETECTOR DERIVED
 *
 * Computed parameters for the detect function
 */

//offsets added to index to move in a specific direction
static int E_O;
static int S_O;
static int W_O;
static int N_O;

//the max/min length of admissible paths
static int LENLIMIT;

/*
 * PROCESSING CONSTANTS
 *
 * Constants for the processing function.
 */

//the number of fields (all integers) that are returned for an ellipse
static const int ELLIPSE_FIELDS = 12;

//TODO initialization not necessary?
static bool THRESHOLDED[MAXSIZE] = { false };

/* DEBUGGING */
#ifdef MOSEY_DEBUG
static char PLOT[MAXSIZE] = { 0 };
#endif


// -------------- UTILITY METHODS ------------- //

// not currently used
//see http://www.azillionmonkeys.com/qed/ulerysqroot.pdf
static unsigned isqrt(unsigned long val) {
    unsigned long temp, g=0, b = 0x8000, bshft = 15;
    do {
        if (val >= (temp = (((g << 1) + b)<<bshft--))) {
           g += b;
           val -= temp;
        }
    } while (b >>= 1);
    return g;
}

//TODO change to a macro?
static int min(int a, int b) {
	return a < b ? a : b;
}

// -------------- DEBUG METHODS --------------- //

#ifdef MOSEY_DEBUG

/* Plots a point with brightness v at x,y */
static void plotPoint(int x, int y, char v) {
  PLOT[y * WIDTH + x] = v;
}

/* Follows length pixels along PATH and plots each with brightness v */
static void plotPath(int x, int y, int length, char v) {
  for (int i = 0; i < length; i++) {
    plotPoint(x,y,v);
    char d = PATH[i];
    switch(d) {
    case E : x++; break;
    case S : y++; break;
    case W : x--; break;
    case N : y--; break;
    }
  }
}

#endif

// -------------- JNI METHODS --------------- //

/* INITIALIZATION */

/* Specifies the width and height of the frame to be processed,
 * along with the degree of intensity averaging to be applied.
 * Calling this function causes all dependent variables to be recomputed.
 * This function may be called any number of times.
 */
static jboolean init(JNIEnv *env, jobject thiz, jint w, jint h, jint o) {

	int width = (int) w;
	int height = (int) h;
	int order = (int) o;
	if (width * height > MAXSIZE) return JNI_FALSE;

	WIDTH = width;
	HEIGHT = height;
	ORDER = order;

	SIZE = WIDTH * HEIGHT;

	HORDER = ORDER;
	VORDER = ORDER;

	HSIZE = 1 << HORDER;
	VSIZE = 1 << VORDER;
	SWIDTH = WIDTH - HSIZE;
	SHEIGHT = HEIGHT - VSIZE;

	HOFFSET = -HSIZE;
	VOFFSET = -VSIZE*WIDTH;

	HHSIZE = HSIZE >> 1;
	HVSIZE = VSIZE >> 1;
	HSWIDTH = WIDTH - HHSIZE;
	HSHEIGHT = HEIGHT - HVSIZE;

	SHIFT = VORDER + HORDER + BIAS_ORDER;

	E_O = 1;
	S_O = WIDTH;
	W_O = -1;
	N_O = -WIDTH;

	LENLIMIT = 4 * (WIDTH + HEIGHT);

	return JNI_TRUE;
}

/* CONVERSION IMPLEMENTATION */

/**
 * Converts values in semi-planar YUV into byted packed RGB.
 */

//TODO make more efficient by tackling pixels in pairs
//see http://blog.tomgibara.com/post/132956174/yuv420-to-rgb565-conversion-in-android
static void convert(jbyte *yuvs, jint w, jint h, jbyte *bytes) {
  //TODO confirm that this cast is safe
  unsigned char *y = (unsigned char*) yuvs;
  unsigned char *off = y + w * h;
  unsigned char *cOff;
  int Y, Cr = 0, Cb = 0;
  int i, j;
  int R, G, B;

  for (j = 0; j < h; j++) {
    cOff = off + (j >> 1) * w;
    for (i = 0; i < w; i++) {
      Y = (int)(*y++);
      if ((i & 0x01) == 0) {
        Cr = (int)(*cOff++) - 128;
        Cb = (int)(*cOff++) - 128;
      }
      R = Y + ((359 * Cr) >> 8);
      if(R < 0) R = 0; else if(R > 255) R = 255;
      G = Y - ((88 * Cb + 183 * Cr) >> 8);
      if(G < 0) G = 0; else if(G > 255) G = 255;
      B = Y + ((454 * Cb) >> 8);
      if(B < 0) B = 0; else if(B > 255) B = 255;
      //TODO assumes little endianess
      *bytes++ = (jbyte) (((G & 0x3c) << 3) | (B >> 3));
      *bytes++ = (jbyte) ((R & 0xf8) | (G >> 5));
    }
  }

}

/**
 * Converts byte-packed RGB values into luminance byte values.
 */

static void luminance(jbyte *bytes, jint len, jbyte *yuvs) {
  //TODO confirm that these casts are safe
  unsigned char *b = (unsigned char*) bytes;
  unsigned char *y = (unsigned char*) yuvs;

  int R, G, B;

  for (; len >= 0; len--) {
    unsigned char j = *b++;
    unsigned char i = *b++;
    //TODO assumes little endianess
    R = i & 0xf8;
    G = ((i << 5) & 0xe0) | ((j >> 3) & 0x1c);
    B = (j << 3) & 0xf8;
    *y++ = (jbyte) ( (R + (G << 1) + B) >> 2 );
  }

}

/* AVERAGING IMPLEMENTATION */

/**
 * Populates the AVERAGE array with luminance values averaged over a square
 * with sides of 2^ORDER pixels. In addition, each pixel is compared against
 * its neighbourhood average and if the pixel is darker, false is recorded
 * in the out array (otherwise true).
 */

//TODO convert implementation to use jint for consistency with JNI
//TODO assumes unsigned ints
void average(char *in, bool *out) {
  //perform horizontal averaging from the original into the averages array
  //this averaging pass could all be done on the original array, but we
  //need to keep that data to then perform the thresholding.
  int index = 0;
  for (int y = 0; y < HEIGHT; y++) {
    int sum = 0;
    for (int x = 0; x < HSIZE; x++) {
      sum += in[index];
      ++index;
    }
    for (int x = HSIZE; x < WIDTH; x++) {
      int z = index + HOFFSET;
      int w = in[z];
      //we divide the sum to get the average later
      AVERAGES[z] = sum;
      sum += in[index] - w;
      ++index;
    }
  }

  //perform vertical averaging from the averaging array back onto itself -
  //we want to preserve the original data (see above).
  index = 0;
  for (int x = 0; x < SWIDTH; x++) {
    int sum = 0;
    for (int y = 0; y < VSIZE; y++) {
      sum += AVERAGES[index];
      index += WIDTH;
    }
    for (int y = VSIZE; y < HEIGHT; y++) {
      int z = index + VOFFSET;
      int w = AVERAGES[z];
      //we make amends for the horizontal summing and bias
      AVERAGES[z] = (BIAS*sum) >> SHIFT;
      sum += AVERAGES[index] - w;
      index += WIDTH;
    }
    index += 1 - WIDTH*HEIGHT;
  }

  //perform the thresholding in 9 separate loops (4 corners, 4 sides, 1 middle)
  //the top HVSIZE rows
  index = 0;
  const int nw = AVERAGES[0];
  const int ne = AVERAGES[SWIDTH-1];
  for (int y = 0; y < HVSIZE; y++) {
    //the first HHSIZE columns
    for (int x = 0; x < HHSIZE; x++) {
      out[index] = in[index] >= nw;
      index++;
    }
    //the middle columns
    const int o = -HHSIZE;
    for (int x = HHSIZE; x < HSWIDTH; x++) {
      out[index] = in[index] >= AVERAGES[x+o];
      index++;
    }
    //the last HHSIZE columns
    for (int x = HSWIDTH; x < WIDTH; x++) {
      out[index] = in[index] >= ne;
      index++;
    }
  }
  //the middle HVSIZE rows
  for (int y = HVSIZE; y < HSHEIGHT; y++) {
    const int yo = (y-HVSIZE)*WIDTH;
    const int w = AVERAGES[yo];
    const int e = AVERAGES[yo+SWIDTH-1];
    //the first HHSIZE columns
    for (int x = 0; x < HHSIZE; x++) {
      out[index] = in[index] >= w;
      index++;
    }
    //the middle columns
    const int o = yo - HHSIZE;
    for (int x = HHSIZE; x < HSWIDTH; x++) {
      out[index] = in[index] >= AVERAGES[x+o];
      index++;
    }
    //the last HHSIZE columns
    for (int x = HSWIDTH; x < WIDTH; x++) {
      out[index] = in[index] >= e;
      index++;
    }
  }
  //the bottom HVSIZE rows
  const int yo = (SHEIGHT-1)*WIDTH;
  const int sw = AVERAGES[yo];
  const int se = AVERAGES[yo+SWIDTH-1];
  for (int y = HSHEIGHT; y < HEIGHT; y++) {
    //the first HHSIZE columns
    for (int x = 0; x < HHSIZE; x++) {
      out[index] = in[index] >= sw;
      index++;
    }
    //the middle columns
    const int o = yo - HHSIZE;
    for (int x = HHSIZE; x < HSWIDTH; x++) {
      out[index] = in[index] >= AVERAGES[x+o];
      index++;
    }
    //the last HHSIZE columns
    for (int x = HSWIDTH; x < WIDTH; x++) {
      out[index] = in[index] >= se;
      index++;
    }
  }

#ifdef MOSEY_DEBUG
for (int i = 0; i < SIZE; i++) {
  PLOT[i] = out[i] ? 100 : 0;
  //graticule
  //if ((i & 0x0f) == 0 || ((i / WIDTH) & 0x0f) == 0) PLOT[i] = 200;
}
#endif
}

/* DETECTION IMPLEMENTATION */

/**
 * Attempts to detect ellipses in the supplied binary thresholded image data
 * by walking contours encountered while walking the pixels from the bottom
 * right-hand corner back. Paths are generated from the contours and then
 * subjected to a number of heuristics before being recorded as an ellipse
 * in the ellipses array. The function returns the number of ellipses detected.
 */

jint detect(bool *data, jint *ellipses, jsize maxEllipseCount) {
  //the number of ellipses we have found
  int ellipseCount = 0;
  //pointer to the next ellipses element to fill
  int *ellipsePtr = ellipses;

  bool *used = USED;
  char *path = PATH;
  int *points = POINTS;

  //set the perimeter to white to sidestep nasty edge problems
  for (int x = 0; x < WIDTH; x++) {
    data[x] = NEUTRAL;
    data[(SIZE-1)-x] = NEUTRAL;
  }
  for (int y = 0; y < SIZE; y+=WIDTH) {
    data[y] = NEUTRAL;
    data[(SIZE-1)-y] = NEUTRAL;
  }

//#ifdef MOSEY_DEBUG
//memset(PLOT, 0, SIZE);
//#endif

  //clear the norths array
  //memset(used, false, SIZE);
  for (int xx = 0; xx < SIZE; xx ++) {
	  used[xx] = false;
  }

  //find a starting point: bottom right hand corner
  int index = SIZE-WIDTH-1;
  //walk backwards so that index matches path coord indices
  for (int y = HEIGHT-2; y >0; y--) {
    index--;
    bool old = NEUTRAL;
    for (int x = WIDTH-3; x > 0; x--) {
      const bool v = data[index];
      //hope that the compiler can remove redundant comparisons
      if (old != NEUTRAL && v==NEUTRAL && !used[index]) {
        int length = 0;
        int i = index;
        char d = BAD_PLACE;
        //assume the path is valid until we show otherwise
        bool valid = true;
        main: do {
#ifdef MOSEY_DEBUG
PLOT[i] = 64;
#endif
          //we calculate an index for the neighbourhood pattern here
          //hope that the compiler can remove redundant comparisons
          int sum = 0;
          if (data[i] == NEUTRAL) sum |= 1;
          if (data[i+1] == NEUTRAL) sum |= 2;
          if (data[i+WIDTH] == NEUTRAL) sum |= 4;
          if (data[i+WIDTH+1] == NEUTRAL) sum |= 8;

          switch (sum) {
          case  1: used[i] = true; d = N; i += N_O; break;
          case  2: d = E; i += E_O; break;
          case  3: d = E; i += E_O; break;
          case  4: d = W; i += W_O; break;
          case  5: used[i] = true; d = N; i += N_O; break;
          case  6:
            if (d == BAD_PLACE) {//we started at a bad place
              valid = false;
              goto skip;
            } else if (d == N) { //TODO revisit this
              d = W;
              i += W_O;
            } else {
              d = E;
              i += E_O;
            }
          break;
          case  7: d = E; i += E_O; break;
          case  8: d = S; i += S_O; break;
          case  9:
            if (d == BAD_PLACE) {//we started at a bad place
              valid = false;
              goto skip;
            } else if (d == E) { //TODO revisit this
              used[i] = true;
              d = N;
              i += N_O;
            } else {
              d = S;
              i += S_O;
            }
          break;
          case 10: d = S; i += S_O; break;
          case 11: d = S; i += S_O; break;
          case 12: d = W; i += W_O; break;
          case 13: used[i] = true; d = N; i += N_O; break;
          case 14: d = W; i += W_O; break;
          default:
        	  LOGE("Invalid local pattern");
        	  return -1;
          }

          path[length++] = d;
          //rare - deals w/ poss. that path is too long
          //we want to continue recording it though so that we don't start walking it again
          if (length == LENLIMIT) {
            valid = false;
            length = 0;
          }
        } while (i != index);
        skip: //we go here when we can't continue with the path
        if (valid) {
          if (length > MINLEN) {
            int count = length >> PTSPACE_ORDER;
            if (count > MAXPOINTS) {
              count = MAXPOINTS;
            } else if (count < MINPOINTS) {
              count = MINPOINTS;
            }
            const int twoCount = count*2;
            int ptIndex = 0;
            int prevJ = -1;
            int nextJ = 0;
            //estimate centroid and accumulate points
            //these record the coordinates of our current location
            int px = x+1;
            int py = y;
            //these are the sums of the x and y coordinates in the path
            int sx = 0;
            int sy = 0;
            for (int j = 0; j < length; j++) {
              if (j == nextJ) {
#ifdef MOSEY_DEBUG
//TODO assumes chars are unsigned
plotPoint(px, py, 255);
#endif
                //compare distances between points - too far = too jagged = not an arc
                if (ptIndex > 0) {
                  int dist1 = abs(px-points[ptIndex-2]) + abs(py-points[ptIndex-1]);
                  int dist2 = j - prevJ;
                  //TODO need to base this on expected curvature?
                  if (dist2 > (dist1*2)) {
                    valid = false;
                    break;
                  }
                }

                //update point data
                points[ptIndex] = px;
                points[ptIndex+1] = py;
                ptIndex += 2;
                prevJ = nextJ;
                nextJ = length*ptIndex/twoCount;

              }
              sx += px;
              sy += py;
              const char pd = path[j];
              switch(pd) {
              case E : px++; break;
              case S : py++; break;
              case W : px--; break;
              case N : py--; break;
              }
            }
            if (valid) {
              //compute the centroid to the nearest pixel
              int cx = (sx+(length>>1)) / length;
              int cy = (sy+(length>>1)) / length;
              //TODO need to deal with zero expectation?
              //check that the shape is convex
              //we also take this opportunity to compute the minimum and maximum distances from the centroid
              //first identify direction of curvature from 3 points
              //we also compute the minimum and maximum distances here
              int x1 = points[twoCount-2];
              int y1 = points[twoCount-1];
              int x2 = points[0];
              int y2 = points[1];
              int x3 = points[2];
              int y3 = points[3];
              //we use the dot product to determine the direction of curvature at each point
              //expected is negative for black and positive for white shapes
              int expected = ((x3-x2)*(y2-y1) + (y3-y2)*(x1-x2)) / NONSIG;
              //TODO could do this in loop that composes point array
              //step over remaining points to confirm consistent direction of curvature (and thus convexity)
              for (int j = 4; j < twoCount; j+=2) {
                x1 = x2;
                y1 = y2;
                x2 = x3;
                y2 = y3;
                x3 = points[j];
                y3 = points[j+1];
                int dot = ((x3-x2)*(y2-y1) + (y3-y2)*(x1-x2)) / NONSIG;
                if (dot != 0) {
                  if (expected == 0) {
                    expected = dot;
                  } else {
                    if (expected < 0 != dot < 0) {
                      valid = false;
                      break;
                    }
                  }
                }
              }
              //check final point triple
              x1 = x2;
              y1 = y2;
              x2 = x3;
              y2 = y3;
              x3 = points[0];
              y3 = points[1];
              int dot = ((x3-x2)*(y2-y1) + (y3-y2)*(x1-x2)) / NONSIG;
              if (expected != 0 && dot != 0 && expected < 0 != dot < 0) {
                valid = false;
              }
              if (valid) {
                //walk the whole path again :-( to
                //compute accurate max & min distance pts
                int minD = INT_MAX;
                int maxD = INT_MIN;
                int mx = -1;
                int my = -1;
                int nx = -1;
                int ny = -1;
                for (int j = 0; j < length; j++) {
                  //compute distance of pt from centroid
                  int d2 = (px-cx)*(px-cx)+(py-cy)*(py-cy);
                  if (d2 < minD) {
                    minD = d2;
                    nx = px;
                    ny = py;
                  }
                  if (d2 > maxD) {
                    maxD = d2;
                    mx = px;
                    my = py;
                  }
                  const char pd = path[j];
                  switch(pd) {
                  case E : px++; break;
                  case S : py++; break;
                  case W : px--; break;
                  case N : py--; break;
                  }
                }

                //TODO perhaps record orientation
                *ellipsePtr++ = cx;
                *ellipsePtr++ = cy;
                *ellipsePtr++ = expected < 0 ? 1 : 0;
                // we report the squares of the distances
                // we avoid the work of calculating their square roots
                // until we are sure we need them
                //*ellipsePtr++ = isqrt(minD);
                //*ellipsePtr++ = isqrt(maxD);
                *ellipsePtr++ = minD;
                *ellipsePtr++ = maxD;
                *ellipsePtr++ = mx;
                *ellipsePtr++ = my;
                *ellipsePtr++ = nx;
                *ellipsePtr++ = ny;
                *ellipsePtr++ = sx;
                *ellipsePtr++ = sy;
                *ellipsePtr++ = length;
                ellipseCount++;
#ifdef MOSEY_DEBUG
//TODO assumes chars are unsigned
plotPoint(cx, cy, 255);
plotPath(x+1,y, length, 127);
#endif
                //EARLY RETURN - we don't have space to record any more info
                if (ellipseCount == maxEllipseCount) return ellipseCount;
              }
            //System.out.println(Arrays.toString(points, 0, ptIndex));
            }
          }
        }
      }
      old = v;
      index--;
    }
    index-=2;
  }
  return ellipseCount;
}

/* PROCESSING IMPLEMENTATION */

/**
 * Processes the frame data (which may be supplied in either yuvsArray as YUV semiplanar data or in bufferObj as RGB_565 data.
 * If RGB data is supplied, the luminance is extracted and written back via the yuvArray.
 * If YUV data is supplied, it is converted to RGB_565 format and stored in bufferObj.
 * The type of frame data is signalled by convertToBuffer: true indicates YUV was supplied, false RGB.
 * After image data has been converted, ellipses are detected, recorded in ellipseArray and their number returned.
 */

static jint process(JNIEnv *env, jobject thiz, jboolean convertToBuffer, jbyteArray yuvsArray, jint w, jint h, jobject bufferObj, jintArray ellipseArray) {
  jbyte *yuvs = env->GetByteArrayElements(yuvsArray, 0);
  jbyte *bytes = (jbyte*) env->GetDirectBufferAddress(bufferObj);
  jsize maxEllipseCount = env->GetArrayLength(ellipseArray) / ELLIPSE_FIELDS;
  jint ellipseCount;

  if (convertToBuffer) {
    convert(yuvs, w, h, bytes);
  } else {
    luminance(bytes, w * h, yuvs);
  }

  if (maxEllipseCount == 0) { //avoid unecessary work
    ellipseCount = 0;
  } else {
    //TODO could simply allocate our own array
    jint  *ellipses = (jint*) env->GetIntArrayElements(ellipseArray, 0);
    average((char*) yuvs, THRESHOLDED);
    ellipseCount = detect(THRESHOLDED, ellipses, maxEllipseCount);
    env->ReleaseIntArrayElements(ellipseArray, ellipses, 0);
  }

  if (convertToBuffer) {
    env->ReleaseByteArrayElements(yuvsArray, yuvs, JNI_ABORT);
  } else {
    env->ReleaseByteArrayElements(yuvsArray, yuvs, 0);
  }

  return ellipseCount;
}

/* SAMPLING IMPLEMENTATION */

/**
 * Samples the thresholded data of batch of x,y coordinates packed sequentially into an pointsArray.
 * The coordinates may safely exceed the bounds of the frame. Sample values are recorded in samplesArray.
 */

static void sample(JNIEnv *env, jobject thiz, jintArray pointsArray, jbooleanArray samplesArray, jint viewWidth, jint viewHeight) {
  if (pointsArray == NULL) return;
  if (samplesArray == NULL) return;
  const jsize pointsLength = env->GetArrayLength(pointsArray);
  const jsize samplesLength = env->GetArrayLength(samplesArray);
  const int length = min( (int)pointsLength / 2, (int)samplesLength);
  jint *points = env->GetIntArrayElements(pointsArray, 0);
  //TODO could simply allocate our own array
  jboolean *samples = env->GetBooleanArrayElements(samplesArray, 0);

  jint *point = points;
  for (int i = 0; i < length; i++) {
    jint x = (*point++) * WIDTH / viewWidth;
    jint y = (*point++) * HEIGHT / viewHeight;
    //TODO is an implicit cast from bool to jboolean safe?
    bool b = x < 0 || x >= WIDTH || y < 0 || y >= HEIGHT ? NEUTRAL : THRESHOLDED[y * WIDTH + x];
    samples[i] = b ? JNI_TRUE : JNI_FALSE;
  }

  env->ReleaseIntArrayElements(pointsArray, points, JNI_ABORT);
  env->ReleaseBooleanArrayElements(samplesArray, samples, 0);
}

/**
 * Samples the average luminance at a batch of x,y coordinates packed sequentially into an pointsArray.
 * The coordinates may safely exceed the bounds of the frame. Sample values are recorded in averagesArray.
 */

static void sampleAverages(JNIEnv *env, jobject thiz, jintArray pointsArray, jintArray averagesArray, jint viewWidth, jint viewHeight) {
  if (pointsArray == NULL) return;
  if (averagesArray == NULL) return;
  const jsize pointsLength = env->GetArrayLength(pointsArray);
  const jsize averagesLength = env->GetArrayLength(averagesArray);
  const int length = min( (int)pointsLength / 2, (int)averagesLength);
  jint *points = env->GetIntArrayElements(pointsArray, 0);
  //TODO could simply allocate our own array
  jint *averages = env->GetIntArrayElements(averagesArray, 0);

  jint *point = points;
  for (int i = 0; i < length; i++) {
    jint x = (*point++) * WIDTH / viewWidth;
    jint y = (*point++) * HEIGHT / viewHeight;
    if (x < HSIZE) x = HSIZE; else if (x >= HSWIDTH) x = HSWIDTH-1;
    if (y < VSIZE) y = VSIZE; else if (y >= HSHEIGHT) y = HSHEIGHT-1;
    //TODO implicit cast to jint
    averages[i] = AVERAGES[(y-HVSIZE) * WIDTH + (x-HHSIZE)];
  }

  env->ReleaseIntArrayElements(pointsArray, points, JNI_ABORT);
  env->ReleaseIntArrayElements(averagesArray, averages, 0);
}

/* DEBUG IMPLEMENTATION */

static jboolean plot(JNIEnv *env, jobject thiz, jobject bufferObj) {
  if (bufferObj == NULL) return JNI_FALSE;
  //TODO should check length of bufferObj
#ifdef MOSEY_DEBUG
  //TODO assumes size of jbyte == size of char
  //TODO assumes little endianess
  //TODO assumes unsigned char
  //convert grayscale plot into rgb565
  char *bytes = (char*) env->GetDirectBufferAddress(bufferObj);
  char *plot = PLOT;
  for (int i = 0; i < SIZE; i++) {
	int v = *plot++;
    *bytes++ = (char) (((v & 0x3c) << 3) | (v >> 3));
    *bytes++ = (char) ((v & 0xf8) | (v >> 5));
  }
  //memcpy(bytes, PLOT, SIZE);
  return JNI_TRUE;
#else
  return JNI_FALSE;
#endif
}

// -------------- JNI REGISTRATION --------------- //

static const char *classPathName = "com/tomgibara/mosey/vision/nat1ve/Native";

static JNINativeMethod methods[] = {
  {"init", "(III)Z", (void*)init },
  {"process", "(Z[BIILjava/nio/ByteBuffer;[I)I", (void*)process },
  {"sample", "([I[ZII)V", (void*)sample },
  {"sampleAverages", "([I[III)V", (void*)sampleAverages },
  {"plot", "(Ljava/nio/ByteBuffer;)Z", (void*)plot },
};

/*
 * Register several native methods for one class.
 */
static int registerNativeMethods(JNIEnv* env, const char* className,
    JNINativeMethod* gMethods, int numMethods)
{
    jclass clazz;

    clazz = env->FindClass(className);
    if (clazz == NULL) {
        LOGE("Native registration unable to find class '%s'", className);
        return JNI_FALSE;
    }
    if (env->RegisterNatives(clazz, gMethods, numMethods) < 0) {
        LOGE("RegisterNatives failed for '%s'", className);
        return JNI_FALSE;
    }

    return JNI_TRUE;
}

/*
 * Register native methods for all classes we know about.
 *
 * returns JNI_TRUE on success.
 */
static int registerNatives(JNIEnv* env)
{
  if (!registerNativeMethods(env, classPathName,
                 methods, sizeof(methods) / sizeof(methods[0]))) {
    return JNI_FALSE;
  }

  //init(176, 144, 3);
  return JNI_TRUE;
}


/*
 * This is called by the VM when the shared library is first loaded.
 */

typedef union {
    JNIEnv* env;
    void* venv;
} UnionJNIEnvToVoid;

jint JNI_OnLoad(JavaVM* vm, void* reserved)
{
    UnionJNIEnvToVoid uenv;
    uenv.venv = NULL;
    jint result = -1;
    JNIEnv* env = NULL;

    LOGI("JNI_OnLoad");

    if (vm->GetEnv(&uenv.venv, JNI_VERSION_1_4) != JNI_OK) {
        LOGE("ERROR: GetEnv failed");
        goto bail;
    }
    env = uenv.env;

    if (registerNatives(env) != JNI_TRUE) {
        LOGE("ERROR: registerNatives failed");
        goto bail;
    }

    result = JNI_VERSION_1_4;

bail:
    return result;
}
